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Abstract* 

In  this  paper  we  pose  a  sequence  of  linear  prediction  problems  that  are  a  little  different 
from  those  previously  posed.  By  solving  the  sequence  of  problems  we  are  able  to  QR  factor 
data  matrices  of  the  type  usually  associated  with  correlation,  pre  and  post-windowed,  and 
covariance  methods  of  linear  prediction.  Our  solutions  cover  the  forward,  backward  and 
forward-backward  problems.  The  QR  factor  orthogonalizes  the  data  matrix  and  solves  the 
problem  of  Cholesky  factoring  the  experimental  correlation  matrix  and  its  inverse.  This 
means  we  may  use  generalized  Levinson  algorithms  to  derive  generalized  QR  algorithms.  /  DTIC 

I  COPY 

which  are  then  used  to  derive  generalized  Schur  algorithms.  All  three  algorithms  are  true  ymsptcTco^ 
lattice  algorithms  that  may  be  implemented  either  on  a  vector  machine  or  on  a  multi-tier 
lattice,  and  all  three  algorithms  generate  generalized  reflection  coefficients  that  may  be  J 
used  for  filtering  or  classification.  ™ 
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I.  Introduction 


In  this  paper  we  pose  a  sequence  of  least  squares  problems  from  the  theory  of  linear 
prediction.  The  problems  are  a  little  different  from  those  originally  posed  in  the  paper  by 
Morf,  et  al.  6  .  The  solution  to  this  sequence  of  problems  produces  QR  factorizations  of  the 
kinds  of  data  matrices  that  are  usually  associated  with  the  covariance,  pre- windowed,  post- 
windowed  and  correlation  methods  of  linear  prediction.  By  QR  factorization  we  mean  the 
computation  of  an  orthogonal  matrix  Q  and  an  upper  triangular  matrix  R  such  that  Y  — 
QR.  That  factorization  is  very  often  used  to  solve  an  overdermined  system  of  equations 
Ya  —  b  in  the  least  squares  sense  by  Ra  -  QTb.  The  results  apply  to  forward,  backward 
or  forward- backward  linear  prediction  problems.  We  interpret  the  QR  factorization  in 
analysis,  synthesis  and  orthogonalizing  terms,  and  use  it  to  generate  Cholesky  factors  of 
the  experimental  covariance  matrix  (  the  Grammian  of  the  data  matrix)  and  its  inverse. 
Then  we  use  the  generalized  Levinson  recursions  derived  by  Friedlander  et  al.  [lj  to 
derive  generalized  recursions  for  computing  the  orthogonal  matrix  in  the  QR  factorization 
of  any  of  the  Toeplitz  (or  concatenation  of  Toeplitz)  matrices  that  can  arise  in  linear 
prediction.  These  recursions  generalize  those  first  discovered  by  Cybenko  1 13]  for  the 
correlation  method  of  linear  prediction.  We  then  use  these  recursions  to  derive  generalized 
Schur  recursions  for  Cholesky  factoring  any  of  the  close-to-Toeplitz  covariance  matrices 
that  can  arise  in  linear  prediction. 

All  of  our  generalized  recursions  are  true  lattice  recursions  that  generate  generalized 
reflection  coefficients.  All  three  algorithms  may  be  implemented  on  a  vector  machine 
or  on  a  multi-tier  lattice.  In  a  procedure  similar  to  that  of  Robinson  and  Treitel  7]  in 
the  scalar  Toeplitz  case,  and  Friedlander  [8  in  the  multi  variable  Toeplitz  case,  we  use 
the  autoregressive  lattice  filter  associated  with  the  generalized  reflection  coefficients  to 
complete  the  QR  factorization  with  the  computation  of  the  upper  triangular  matrix. 

Finally,  w’e  show  how  our  generalized  Levinson,  QR,  and  Schur  algorithms  may  be 
extended  to  multivariable  (or  vector)  linear  prediction  problems.  The  multivariable  results 
may  be  applied  to  a  variety  of  multidimensional  problems,  as  wrell. 

II.  Least  Squares  Problems  in  Linear  Prediction 

Let  y  =  [y0,  3/i ,  •  •  • ,  y/v-i  \T  denote  an  N  sample  snapshot  of  the  stationary  time  series 


{yt}.  From  this  snapshot,  we  would  like  to  identify  an  autoregressive  or  whitening  model 
for  the  time  series  {y(}.  This  model  takes  the  form 


«  =  K  ;  a"  -  l- 


(II- 1) 


where  {u”}  is  a  white  sequence  with  zero  mean  and  variance  The  interpretation  is 
that  the  digital  filter  J4n(z)  =  V]”=0a"z  '  whitens  the  time  series  {y,}.  The  whitening 
model  may  be  written  as  the  predictor  model 


Vt  =  yiJt 


(II  2) 


The  variance  of  the  white  sequence  {u<}.  or  equivalently  the  mean  squared  error  between 
yt  and  the  one-step  ahead  predictor  yt ,  is  a‘n. 

Prediction  :  Our  procedure  for  identifying  a  model  „4.n(z)  will  be  to  form  a  sequence 
of  predictions  of  the  form 


=  :  a*  =  1 


(II  3) 


and  to  let  the  predictor  order  range  from  p  -  0  to  p  =  n.  The  error  between  y,  and  the 
p order  predictor  yf  is,  of  course, 


uf  =  yt  -  =  5^  afyi-. 


VI*) 


We  shall  be  interested  in  a  window  of  these  errors  for  which  the  time  index  i  satisfies  the 


condition 


k  <  i  +  (n  -  p)  <  l. 


(II- 5) 


As  the  predictor  order  increases  from  p  —  0  to  p  —  n,  the  window  of  length  (/  -  k  +  1) 
moves  from  left  to  right  across  the  data  set,  as  illustrated  in  Figure  1.  The  indexes  k  and 
l  may  be  chosen  to  select  among  the  various  techniques  of  linear  prediction.  The  data 
values  outside  the  range  10, . . . ,  JV  —  1]  are  set  to  zero.  In  the  covariance  method  of  linear 
prediction,  k  =  n  and  l  —  N  ~  1,  in  the  correlation  method,  k  ~  0  and  /  =  N  —  1  +  n,  as 
illustrated  in  Figure  2. 


Let  us  write  out  the  error  equations,  over  the  window  just  defined,  for  the  pin  order 
predictor  as  follows  : 


yi  i  yi 


ll  n  -  p 


1  f  n  -  p 


(II- 6) 


The  compact  notation  is 


Vk  4  1 

Vki  Vk 


Y  -  =  r  4P  =  up 

0 


1  fc -  n -  p-K  1 

,p 

k  —  n  —  p 


This  scheme  may  be  reproduced  for  p  -  0  to  p  =  n  to  obtain  the  set  of  equations 


5.4  =  V 


(II- ?) 


(II- 8) 


where  5'  is  the  Toeplitz  data  matrix  just  defined,  and  the  matrices  .4  and  U  are  given  by 


=  U",...,^lp, . 4n 


(II.  9) 


/  -  n  -  p 


T:",...,  rp . rn 


(77.10) 


uk-n^  1 

,,o 


1  Jfc  —  p 


Least  Squares  :  The  Grammian  of  the  error  matrix  V  is 


UTU  =  ATYTY  A 


(77.11) 


■SK 
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The  pp  element  of  the  Grammian  is 


a2p  =  (UP)TUP  -  (ap)r0T  YTY  g  (Ap )TYTYAP, 


(II- 12) 


which  is  just  the  accumulated  squared  prediction  error  over  the  window  defined  for  the 
predictor  of  order  p.  This  squared  error  may  be  minimized  with  respec  t  to  the  coefficients 
ap  under  the  constraint  that  (ap)T b  -  ap  1.  The  appropriate  regression  equation  is 


>  (k°,’)Vttv  -  a(s nTt)  =  <> 


(II-  13) 


The  solution  for  ap  is  then 


!  i  o  yty  * 


-cr2J  ■ 


(II- H) 


the  matrix  \I  Oi  is  just  there  to  reduce  the  length  of  the  vector  to  1.  When  this  solution 
is  written  out  for  p  =  0, 1, . . . ,  n,  the  result  is 


Y7YA  =  LD 2 


(11.15) 


where  L  is  a  lower  triangular  matrix  with  ones  on  its  main  diagonal  and  D 2  is  the  diagonal 
matrix  containing  the  prediction  errors  : 


D2  -  Diag  \<rl<r*,...yn} 


(II- 16) 


QR  Factors  :  The  equation  YTYA  -  LD2  characterizes  the  least  squares  predictors 
for  p  -  0  to  p  =  n.  The  right  hand  side  is  lower  triangular.  If  both  sides  of  the  equation 
are  pre-multiplied  by  AT .  then  the  left  hand  side,  namely  ATYTY A  is  symmetric  and  the 
right  hand  side  is  lower  triangular.  So  the  right  hand  side  must  be  diagonal.  This  means 
the  least  squares  solution  for  the  columns  of  A  produce  the  following  equations  : 


YA  =  V 


AtYtYA  =  VTU  =  D2 


(ii.\r 


It  follows  that  YA  =  U  is  a  QR  factorization  of  the  data  matrix  Y  when  ,4  solves  the 
sequence  of  least  squares  problems  we  have  posed.  This  QR  factorization  (really  a  Gram- 
Schmidt  orthogonalization  of  Y )  may  be  rewritten  as 


Y  =  UHJ 


(II- is) 


! 


where  HT  is  upper  triangular.  Then  the  QR  factorization  may  be  given  the  following 
analysis,  synthesis  and  orthogonalization  interpretations  : 

i)  YA  —  U  :  Y  is  analyzed  or  whitened  by  the  analysis  matrix  .4  to  produce  the  white 
matrix  V  (UTU  -  D2). 

ii)  Y  =  VH 7  :  Y  is  synthesized  from  the  white  matrix  V  by  the  upper  triangular  synthesis 


matrix  HT  -  A  1 . 


iii)  YtU  ~  HD 2  :  the  data  matrix  Y  and  the  white  matrix  U  are  causally  orthogonal, 
with  H  describing  the  cross-correlation  between  the  input  1  and  the  output  V . 
Second-Order  Equations  :  Let  us  define  the  experimental  correlation  matrix  R  to  be 


the  Grammian  of  the  data  matrix  Y  : 


R  =  YtY 


(II- 19) 


Then,  from  the  analysis  and  synthesis  equations  we  can  interpret  the  analysis  matrix  A  as 
a  Choleskv  factor  of  R _1  and  the  synthesis  matrix  H  as  a  Choleskv  factor  of  R  : 


i)  AT R A  =  D2  :  A  is  a  Cholesky  factor  of  R  1 . 

ii)  R  =  HD2  Ht  :  H  is  a  Cholesky  factor  of  R. 

iii)  RA  ~  HD 2  :  .4  analyzes  R  to  produce  H. 


Backward  equations  :  The  prediction  operation  may  also  be  written  backwards  in  time 
(with  respect  to  future  values  of  the  time-series).  The  one-step  backward  model  becomes 


b?yt  +  i  -  : 


K  =  1. 


(77.20) 


where  {r"}  is  a  white  sequence  with  zero  mean  and  variance  t2.  As  before,  lower  order 
predictors  are  introduced  to  form  a  family  of  order  increasing  predictors 


Kyt+l  =  vf  ;  iff  =  1. 


(77.21) 


for  p  —  0  to  p  —  n.  The  prediction  errors  are  written  for  the  time  index  satisfying  the 
following  condition  : 


k  <  t  +  p  <  /. 


(77.22) 


The  error  equations  for  the  order  p  predictor  may  be  written  as 


yk  Pk -  1  ■  ■  •  yk -  n 

- 

''Ip 
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ry 

vl  p-] 
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The  data  matrix  X  in  this  equation  is  related  to  Y  by  the  exchange  matrix  .7 
The  compact  notation  is  then 


(77.23) 


.Y  =  JYJ. 


X 


y 

y 


XBP  -  Vp 


(77.24) 


This  scheme  may  be  reproduced  for  p  —  0  to  p  —  n  to  obtain  the  set  of  equations 


XB  ^  V  (77.25) 

B  -  iR°, . . . ,  Bp, . . .  ,Bn  V=.  (77.26) 

The  least  squares  solution  under  the  constraint  -  1,  when  written  out  for  p  -  0  to 
p  =  n,  leads  to  the  matrix  equation  : 

B‘XtXB  =  VTV  --  E2  --  Diag  [t2.,..  ..r2p . r2  (77.27) 

It  then  follows  that  X  -  YB~'i  —  YGr  is  a  QR  factorization  of  A .  and  the  second  order 
Choleskv  equations  for  the  Grammian  XTX  -  J RJ  are 

BTXTXB  =  E 2  XT X  =  GE2Gt  (77.28) 

The  connections  with  the  forward  case  come  from  the  relationship  between  A  and  i  . 


The  QR  factorization  of  X  is  a  “QL”  factorization  on  Y.  The  upper-lower  Choleskv 
factorization  of  (A'^A*)'1  is  a  lower-upper  Choleskv  factorization  of  (YTY)  1  and  the 
lower-upper  Choleskv  factorization  of  XTX  is  an  upper-lower  Choleskv  factorization  of 


l-l'l  ***'•«•  'ftl  M  *4|  '»u 


i.|  *.»«*»  **4  «.ti 


Forward- Backward  equations  :  When  the  time  series  is  known  to  he  stationary,  its 


statistical  properties  are  not  modified  when  time  is  reversed.  This  property  may  he  very 
important  for  some  applications,  and  the  family  of  models  may  he  forced  to  satisfy  this 
property.  In  other  words,  the  forward  and  backward  prediction  errors  are  written  with 


respect  to  the  same  predictor  coefficients  : 


^  afy,  i 


E°r: 


<  i. 


(77. 29) 


The  compact  notation  is 


Y  qr 
X  0 


A*  -  JYJ 


(77.30) 


When  these  equations  are  written  out  for  prediction  orders  from  p  0  to  p  n.  we  have 
the  coupled  prediction  equations 


(77.31) 


The  Grammian  of  the  error  matrix  is 


j'T  \’T  [ 


T  vT  vT 


A*  Y1  X 


A  --  At\YtY  -  XTX  A 


(77.32) 


The  pp^  element  of  this  Grammian  is 


al  =  (AP)T\YTY  +  XT X  Ap 


(77.33) 


which  is  just  the  accumulated  squared  prediction  errors  in  the  forward  and  backward 
predictions. 


As  before,  the  solution  of  the  family  of  least  squares  problems  produces  a  matrix  .4  that 
orthogonalizes  the  data  matrix  ^  .  and  Cholesky  factors  the  Grammian  Y1  Y  -  Xr  X. 

Linear  Statistical  models  :  There  is  one  more  useful  comment  to  be  made  about  the 


equations  of  linear  prediction.  Suppose  y'  is  the  snapshot  of  a  data  set.  It  might  be 
obtained  from  a  multi-sensor  array  or  from  overlapped  windows  of  a  time  series.  This 
snapshot  may  be  analyzed  into  a  white  vector  u1  by  fitting  a  linear  prediction  or  analysis 


model  : 


ATy'  -  u* 


(77.34) 


where  Ar  is  a  lower  triangular  matrix  with  unit  diagonal  elements.  Call  R  the  experimental 
correlation  matrix  and  D 2  the  experimental  analysis  matrix  : 


Vy*  (y')T  R  Yu'fu’)7  D~ 


(11. 35) 


Then  the  connection  between  J?  and  D 2  is 


ArRA =  D- 


(11. 36) 


If  the  columns  of  .4  are  chosen  to  minimize  the  diagonal  elements  of  D‘.  term  by  term, 
then  D 1  is  diagonal,  and  the  factorization  obtained  is  a  Cholesky  factorization.  Similarly 
the  synthesis  model  mav  be  written 


y'  r-  Hu'  Ht  -  A 


(II  .37) 


The  connection  with  the  prediction  equations  comes  with  the  choice  of  y1  as  the  column 
of  YT .  the  Toeplitz  data  matrix.  With  such  a  choice  for  y\  the  linear  statistical  model  is 
really  another  way  to  write  down  the  whole  set  of  increasing  order  prediction  problems. 
When  the  matrix  R  is  replaced  by  the  expectation  of  y'(y')T .  then  this  formalism  reduces 
to  that  of  3  . 

Summary  :  By  solving  the  right  sequence  of  least  squares  problems,  we  QR  factor 
a  data  matrix  and  produce  Cholesky  factors  of  the  experimental  covariance  matrix  and 
its  inverse.  This  is  fundamental.  We  can  either  think  of  a  QR  factoring  of  the  data 
matrix  as  a  square  root  method  of  factoring  the  experimental  correlation  matrix  and  its 
inverse,  or  we  can  think  of  a  Cholesky  factoring  of  the  experimental  covariance  matrix 
and  its  inverse  as  a  square  method  of  obtaining  the  “R"  part  of  a  QR  factorization.  We 
can  also  think  of  the  QR  factor  T.4  -  l\  the  Cholesky  factor  A1  RA  D2 .  and  the 
Cholesky  factor  R  HD2HT  as  three  different  ways  of  characterizing  the  matrix  .4  (or 
its  inverse)  that  contains  order-increasing  prediction  filters.  In  these  characterizations,  the 
Toeplitz  structure  of  T,  and  the  close-to- Toeplitz  structure  of  R.  may  be  used  to  derive 
fast  algorithms  for  computing  A  (or  its  inverse).  In  these  algorithms,  generalized  reflection 
coefficients  are  used  to  update  recursions. 


In  section  III  of  this  paper,  we  review  how  the  generalized  Levinson  recursions  derived 
by  Friedlander  et  al.  I  produce  a  fast  algorithm  for  computing  A  from  the  factorization 
A1  It  A  D2 .  These  recursions  are  used  in  section  I\  to  produce  a  fast  algorithm  for 
computing  the  orthogonal  matrix  l  from  the  factorization  }’.4  {'.  Then  the  recursions 

for  l  are  used  in  section  V  to  produce  a  fast  algorithm  for  H  from  the  factorization 
R  H D2  H1 .  All  these  algorithms  are  true  lattice  algorithms  that  produce  generalized 
reflection  coefficients.  In  this  way  we  will  have  derived  fast  QR  algorithms  for  any  of  the 
data  matrices  that  arise  in  linear  prediction  problems  and  generalized  Schur  algorithms 
for  any  of  associated  experimental  covariance  matrices. 

III.  Factoring  7?  1  into  its  Cholesky  factors 

The  problem  of  factoring  R  1  is  the  problem  of  finding  .4  in  the  diagonalization  : 


AT R A  D 2  Diag  a2 .  er2 


(III  A 


This  equation  may  be  written  as 


and  read  out  as  follows  : 


7?.4  HD 2  ;  H  -  A 


R,a'  a2fit  a2  0. 


(Ill  .2) 


(111.2) 


\Y  here  R,  is  the  (t  •  1)  by  (it  1)  top  left  submatrix  of  R.  When  t  is  incremented  to  7  ■  1 
then,  of  course,  a  new  column  and  a  new  row  are  added  to  7?,.  If  the  resulting  matrix 
7?,.  ]  has  a  simple  recursive  dependence  on  77,,  then  there  is  reason  to  hope  for  a  recursive 
dependence  of  a '  1  on  a'.  This  was  the  insight  of  Friedlander  et  al.  1  . 

The  matrix  7?  is  YJ }  (or  A’7  A',  or  )  1  }  -  A  1  A  )  with  )  and  A'  Toeplitz.  This  means 
7?  is  close-to-Toeplitz.  Let  us  denote  the  (ti  -*  1)  by  (ti  -  1)  symmetric,  non  negative  definite 
correlation  matrix  7?  as  follows  : 


S 


7u.ll  ro.i 

7].u  r,j 


1  ,U 

rn  \ 


7(1. n  1  n 


7n  1  ,n 

7n  ,n  1  7n  r, 


(777.4) 


The  shifted  difference  matrix  T  is  the  n  bv  n  matrix  b  R  : 


1 


{III. 5) 


The  rank  of  b  R  is  the  displacement  rank  a.  We  stress  the  cases  where  0  <  o  <  4.  as  these 
correspond  to  the  least  squares  problems  of  linear  prediction.  The  decomposition  of  b  R 


mav  be  written  1 


b R  = CZCT 


(III. 6) 


where  C  is  an  n  by  a  matrix  and  E  is  an  a  by  a  diagonal  signature  matrix,  containing  -  1 
or  -1  on  its  diagonal. 

The  fundamental  equation  used  in  the  derivation  of  fast  algorithms  is  the  update  for 
the  matrix  Rt,  using  Ct.  which  consists  of  the  first  (i  -  1)  rows  of  C  : 


R,*i  - 


ro ,»-» i 


rt-»l,0  ■■■ 

r»M  !rU,l  '  '  '  r0,i4  1 

fj,0  i 


f,-n,0  1 


{III. 7) 


c,sc,7 


The  idea  here  is  to  correct  a  Toeplitz  approximation  for  R, .  j  with  a  low  rank  matrix 
C,EC7.  Note  that  R  -  Rn  and  R y  =  ru  y. 

When  R  is  Toeplitz  then  r,^  -  r  t  }  .  which  means  that  b  R  is  zero.  In  the  more 
general  case,  5  7Y  has  a  displacement  rank  equal  to  zero  in  the  correlation  method  of 
linear  prediction,  one  in  the  pre-  and  post-windowed  methods  of  linear  prediction,  two 
in  the  covariance  case  of  linear  prediction  and  four  in  the  forward-backward  covariance 
method  of  linear  prediction.  Intermediate  forward-backward  methods  with  a  displacement 
rank  of  two  may  be  introduced,  by  using  forward-backward  methods  in  either  the  pre 
windowed  or  the  post-windowed  methods  of  linear  prediction. 

Let  a ’  -  la’,...,a’,l  7,  denote  the  first  (?  -+-  1)  elements  of  .4'.  the  column  of  .4. 

From  the  Choleskv  decomposition  R~]  -  AD  2AT.  it  follows  that  RA  HD2,  where 


i«!i 


A  T .  Read  the  column  of  this  equation  : 


If  a1'1  is  approximated  bv  .  then  the  (?  -+■  1)  version  of  this  equation  is 

a* 


(7/7.8) 
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which  mav  be  written  as 


R  ° 
R"'  a' 


•  <  -  •  />, 
0  :  C, 

1  0 


(111.10) 
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Define  the  (i  -*  1)  by  (o  -+  1)  matrix  a'  as 


1  0' 

0 

:  C,., 

o 


(777.12) 


where  A,  is  an  (o  *  1)  by  (o  -  l)  error  matrix,  so  that 
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(777.14) 


SSsSSskk 


raw 


The  prediction  error  cr* ,  and  the  error  matrix  A,  are  updated  as  follows  : 

*1,  -  ^2(1  Kh  A,  '  A'.-i) 

trf,  ,  A,T  1  -  ct*(A,  Am  1  A',7  ,  ) 


(7/1.20) 


Note  that  these  two  equations  are  the  same  if  A,  1  and  A’,_  i  is  a  scalar,  and  that  A, 
never  needs  to  be  inverted,  as  only  A,  1 K,^  \  is  used  in  the  recursions. 

IV.  Factoring  )'  into  its  QR  factors 

Using  the  recursions  for  the  columns  of  A,  we  propose  to  find  the  corresponding  re¬ 
cursions  for  the  columns  of  the  orthogonal  matrix  U,  using  the  QR  equation  }\4  =  V.  For 
the  correlation  method  of  linear  prediction  we  use  the  technique  of  Rialan  and  Scharf  ;4 . . 
For  all  the  other  cases  of  linear  prediction  we  generalize  their  technique. 

Correlation  method  of  LP  :  In  the  correlation  method  of  linear  prediction,  k  -  0  and 
1  -  N  1  -  n.  which  means  that  the  data  matrix  Y  looks  like  this  : 


2/A  1  •  •  •  2/ A'  -  n  2/A-n  ]  •  •  •  3/0  0 


3/0  0 

3/i  3/o 


(jr.i) 


Then  the  correlation  matrix  YTY  =  R  is  Toeplitz,  symmetric,  and  positive  semi-definite. 

The  Cvbenko  recursions  10'  on  the  columns  of  V  may  be  derived  by  using  the  Levinson 
recursions  on  the  columns  of  A  to  induce  recursions  on  the  columns  of  L  ,  using  1.4  V . 
Define  Y,  to  be  the  first  (i  -  1)  columns  of  Y  and  U'  to  be  the  column  of  V.  Then 


introduce  an  auxiliarv  vector  V  .  so  that 


i\4'  -  v.q'  =  u*  v.4  =  y,s'  v 


(1Y.2) 


We  saw  in  section  II  that  the  column  of  V  contains  the  forward  prediction  errors  of 
order  i.  Similarly  V  contains  the  backward  prediction  errors  of  the  same  order.  Reproduce 
these  equations  for  (i  +  1)  and  use  the  Levinson  recursions  to  get 


y.41^1  =  v" 1  -  y  (z.4’  -  *•,. 


(IV-  3) 


As  we  have  l'Z.4'  -  ZYA'  (as  long  as  the  last  element  in  .4'  is  equal  to  zero),  then 


t/,4)  ZU'  ■+  k,^U 


(IV-  4) 


Doing  for  the  same  for  V  leads  to  the  following  recursions  : 

U -  V  -  k,. ,  2  U‘ 

V,  }  Z  V'  -  fr„,  r'  for  /  0 . N  1. 


(/r.5) 


These  recursions  are  initialized  using 


T"  -  1’  -  y\-  i  ...  Pi  Pu  0 


*2  -  (tr0)rr(l 


(/r.6) 


]  is  the  reflection  coefficient,  and  is  computed  from  the  internal  variables  by  the  follow 


ing  inner  product 


k,^  a;  =  Z  U' 


This  comes  from  the  familiar  computation  of  the  reflection  coefficient 

,  a]  -  -  j  ro  C]  ...  rn  1  Z  A' 

=  -  ( y\  - 1  •••  Pt>  o  •  •  •  0  Y  Z  A' 

=  -(r°)7  z  iTi 


(iv.  s) 


Alternates  formulas  for  the  reflection  coefficient  are  given  by  the  fact  that  l7,+  1  is  orthog¬ 
onal  to  the  previous  V1  (j  =  0 as  well  as  the  previous  V 1  : 


(r’)Tzr*  (rYzr 


(tM)rr 


[V)TV 


Note  that  <r2  -  (l  ')T  l  '  =  (['  )Tf  .  This  algorithm  for  computing  the  orthogonal  matrix 
in  the  QR  factorization,  is  composed  of  parallel  updates  for  the  vectors  f  and  V  ,  and  a 
side  computation  containing  an  inner  product  for  fr,,i<r2,  and  the  update  for  cr2_,  : 


<£i  -  "?(1  - 


(IV.9) 


The  Cybenko  algorithm  may  be  extended  to  more  general  matrices  1  .  corresponding 
to  the  other  cases  of  linear  prediction,  using  the  generalized  Levinson  recursions.  The 
only  difference  lies  in  the  dimension  (a  4  1)  of  the  auxiliar  matrix  A  .  which  induces  a 
different  dimension  in  the  matrix  V  ■  We  saw  that  the  very  special  shape  of  V  in  the 
correlation  method  of  linear  prediction  allowed  us  to  permute  factors  in  Y Z A' .  The  loss 


B 


of  the  upper  triangle  of  zeros  in  }’  in  the  covariance  method  forbids  us  from  making  the 
same  manipulation.  We  therefore  develop  a  “generic"  algorithm  based  on  the  extension 
of  1  with  the  upper  triangle  of  zeros  (when  it  is  missing).  The  matrix  V  in  the  QR 
factorization  of  V  for  the  covariance  method  of  linear  prediction  will  then  be  contained 
within  the  orthogonal  matrix  of  the  “generic"  case.  The  pre  and  post-windowed  methods 
follow  as  special  cases. 

In  the  correlation  method  of  Linear  Prediction,  the  matrix  H  may  be  computed,  if 
needed,  by  using  the  "impulse”  response  and  internal  variables  of  the  AR  (or  recursive) 
lattice  as  described  in  |3  .  In  section  VI.  we  will  show  how  this  procedure  generalizes  in 
the  other  cases. 

Generic  case  :  We  suppose  now  that  we  have  the  following  Toeplitz  matrix  IT  : 


'  yi 

o  ...  o  - 

yi~  t 

yi  ’ '  ■  : 

IT  = 

'• 

0 

yi  n 

yi 

-  Vk  -  n 

Vk- 

The  data  matrix  Y  has  QR  factor  } \4  - 

V.  Therefore 

T 

Y 


(jr.io) 
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which  means  the  orthogonal  matrix  V  is  embedded  in  a  larger  matrix.  Define  the  i 
column  of  IT. 4  as  follows  : 

\YA' 


W 

V' 


Similarly 


H\4 


U 

—  I 

V 


(IV- 12) 

(7V.13) 


As  in  the  correlation  case,  we  reproduce  these  equations  for  ( )  t  1 )  and  use  the  generalized 
Levinson  recursions  for  .4*  and  .4  to  get 
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(IV- 14) 


16 


The  computation  of  Kl+i  depends  then  on  the  method  used. 

Covariance  method  of  LP  :  For  the  covariance  method  of  linear  prediction, 
I  N  1.  Then  the  shifted  difference  matrix  6  R  has  rank  2  : 

„  T 
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y-n  Vn  1 
Vn  -+■  J  Vn  J 
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Write 


r.4*  =  =  u'  ya  =  r,s*  =  u' 

with  U  an  (N  -  n)  by  3  matrix.  The  “generic”  recursions  may  be  used  for  U' 
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The  initialization  is 
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The  three  components  of  the  reflection  coefficient  itl4]  are  computed  using 


(7V 


which  comes  from  the  use  of 


-  -(ZA*? 


yu  y\-n 


,0  2/ri  1  y.\  1 


(71  .22) 


Only  one  inner  product  is  necessary  per  update.  Note  that  f'(A'  n)  is  the  last  value 
contained  in  the  vector  Lr‘,  and  that  U'( n)  is  the  last  value  of  the  vector  Ul .  Reflection 
coefficients  are  computed  without  first  computing  correlations. 

An  alternate  formula  may  also  be  obtained  as  in  the  correlation  method,  using  the 
orthogonality  property  of  the  columns  of  U  (but  is  much  more  expensive)  : 

»  =  [o  (fY]z["!]-[o  (C')7j  ["■]  V*-. 

At  this  point,  if  one  needs  the  whitening  matrix  .4  (or  its  last  column,  the  coefficients  of 
4„(:)).  then  the  generalized  Levinson  recursions  may  be  used  together  with  the  reflection 
coefficients  computed  by  the  QR  recursions.  In  this  approach,  what  is  really  done  is  to 
compute  the  impulse  response  of  the  lattice  filter  associated  with  these  reflection  coefficients 
(see  section  VI). 

Post-windowed  method  of  LP  :  In  the  post -windowed  method  of  linear  prediction,  k 
n  and  1  -  N  4  n  -  1.  Then  Y  and  the  shifted  difference  matrix  t  It  are  defined  bv 


' y\ ■  - 1  yr\-2  •  •  •  y\~n  i 

0  y*  i 

.0  ...  0  y.v  , 

C  -  y„  ...  y„-  ]  ]T  and  V  - 


(JV.23) 


(IV-  24) 


The  displacement  rank  o  is  unity  which  means  that  the  generalized  Levinson  algorithm 
will  use  an  n  by  2  vector.  The  “generic"  recursions  may  be  simplified  in  this  ease  to  : 


17*41  =  Z  U'  -  U  Ai  1  Kl+^ 
u'+ ’  =  &'  +  Z  U'  K*, 


(IV  25) 


for  i  -  0 . n  1 . 
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The  initialization  is 
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Now  the  reflection  coefficient  A\+i  is  given  by 

*]Klx  !-(t'u)T  z  r\  r*(Ar) 


(JV.27) 


where  Cr,(.V)  is  the  last  element  of  t’1.  Note  that  the  side  computations  involve  an  inner 
product,  and  the  updates  for  A,  and  cr~ .  as  in  the  corresponding  generalized  Levinson 
algorithm. 

Pre- windowed  method  of  LP  :  For  the  pre- windowed  method  of  linear  prediction,  k  - 
0  and  /  .V  -  1.  Then  b  R  is  rank  1.  and  Y  and  b'R  are  defined  bv 
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(7V.28) 


(/r.29) 


Note  that  the  pre-window  matrix  V  is  related  to  the  post-window  matrix  Y  by  the  exchange 
matrix  J .  JYJ  has  the  same  triangle  of  zeros  on  the  top  as  in  the  post-windowed  method. 
The  only  difference  lies  in  the  fact  that  time  is  reversed,  or  in  other  words  y,  is  replaced  by 
y\  i  ,.  The  fast  algorithm  for  U  in  the  post- windowed  case  may  then  be  used  on  JY J . 

Backward  Linear  Prediction  :  In  the  backward  method  of  linear  prediction,  the  only 
difference  with  the  forward  recursions  is  that  the  entries  in  the  data  matrix  V  are  altered. 
The  shape  of  Y  remains  unchanged.  In  other  words,  time  is  reversed.  This  means  that 
the  fast  algorithms  for  Y  may  be  used  for  V  with  just  a  modification  of  the  initialization 
of  the  vectors. 

Forward- Backward  Linear  Prediction  :  In  the  forward-backward  method  of  linear  pre¬ 
diction,  the  data  matrix  is  the  concatenation  of  two  Toeplitz  matrices.  For  the  correlation 
method,  YTY  -  JYTYJ  meaning  that  the  forward-backward  extention  leads  to  the  same 
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Cholesky  factors  .4  and  if,  and  D 2  is  just  multiplied  by  a  factor  of  two.  In  the  covariance 
method,  the  forward- backward  extension  leads  to  the  so-called  modified  covariance  method 
of  linear  prediction,  and  one  way  to  compute  the  highest  order  predictor  -4n(z)  is  to  use 
the  algorithm  of  Marple  13  .  But  the  lower  order  predictors  introduced  in  13  do  not 
solve  the  same  set  of  least  squares  problems  that  we  have  defined.  We  introduce  in  the 
Appendix  an  alternate  to  this  algorithm  to  find  a  characterization  of  the  highest  order 
predictor,  based  on  the  generalized  reflection  coefficients  computed  with  the  recursions  for 
the  columns  of  the  orthogonal  matrix.  The  impulse  response  of  the  associated  lattice  filter 
is  then  equal  to  the  highest  order  predictor. 

V.  Factoring  R  into  its  Cholesky  Factors 

The  LU  factorization  of  R  may  be  written  R  —  HD2HT.  In  the  Toeplitz  case,  the 
LeRoux-Gueguen  algorithm  2|  may  be  used  to  compute  H  directly.  In  all  the  linear 
prediction  cases,  or  in  other  words  when  R  -  YTY ,  the  recursions  for  the  columns  of  H 
are  easily  induced  from  the  QR  factor  recursions  by  simply  premultiplying  the  recursions 
by  YT.  This  was  first  done  in  [4j  for  the  Toeplitz  case. 

Denote  by  H'  the  vector  given  by  YTU *  —  H'.  and  by  H  the  matrix  given  by 

i  ^  i  #  — » t 

Y 1  U  H  .  Use  the  recursions  for  V'  and  V  to  obtain  4 he  recusrions  : 


H  -t YT  Z  U'Kj+i 


ri4  1  _  \-T 


Y‘  Z  V'  +  H  A;1  A, 


Set  the  first  row  of  H  to  zero  to  be  able  to  replace  YT Z V'  by  Z H’ .  The  algorithm  consists 
then  of  the  recursions 

h'"1  -  h’  +  Z  H’Kl , 

(V.2) 

H ,  +  1  -  Z  if*  4  H  A", 4 j 


A  more  general  approach  to  this  derivation  consists  of  using  the  generalized  Levinson 
recursions  directly  to  induce  recursions  for  the  columns  of  H .  Then  there  is  no  need  for  R  to 
equal  YTY .  This  algorithm  is  a  generalization  of  the  vector  version  of  the  Leroux-Gueguen 
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algorithm  derived  in  j3i.  Define  h'  and  g'  by 
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Note  that  h'  is  a  vector  of  length  (n  i  1)  and  g'  is  an  (n  -  i  1)  by  (a  4  1)  matrix. 
Our  goal  is  to  derive  recursions  for  h'  and  g' .  without  carrying  the  update  equations  for 
.4'  and  A  .  Some  simplifications  arise  with  the  definitions  of  H'  and  H 
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The  coupled  recursions  are  : 
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for  j  -  0 . n  1. 
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Note  that  H'  is  the  column  of  HD 2,  or  similarly  H  ‘  2  -  the  column  of  H . 
Note  also  that  the  first  i  elements  of  H'  and  the  first  (  j  -  1 )  rows  of  H  are  equal  to  zero. 
The  coefficient  Kl+ 1  may  be  read  out  of  the  recursions  as 
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so  that  Kj, , a2  equals  the  first  non-zero  row  in  H  .  These  coupled  recursions  include 
the  update  for  <r2  -  H'(i).  the  first  non-zero  element  in  H'.  The  only  necessary  side 
computation  is  the  update  of  A,  : 


j  A, .  ,  (7-,2(  A,  K,,  ]  A',7  j ) 


This  algorithm  generalizes  the  LeRoux-Gueguen  algorithm  to  close- to-Toeplitz  matri¬ 
ces  and  places  the  recursions  in  a  vector  format.  The  algorithm  is  fixed  point  and  no  inner 
product  is  necessary  to  compute  the  reflection  coefficients.  The  recursions  are  equivalent 
to  the  lattice  algorithm  of  Friedlander  5  .  The  only  differences  lie  in  the  normalization  of 
the  lattice  and  in  the  organization  of  the  cells. 

If  the  reflection  coefficients  are  known,  it  is  a  challenging  problem  to  see  if  and  how  the 
correlation  coefficients  and  the  matrix  H  may  be  computed  from  them.  In  the  correlation 
case.  Robinson  and  Treitel  ,7!  solved  this  problem  by  observing  that  the  all-pole  lattice 
filter  has  an  output  equal  to  the  causal  part  of  the  correlation  sequence  when  the  input 
is  zero  and  the  state  is  initialized  at  r„,0 . 0  .  The  multi-channel  case  was  studied  by 


Friedlander  8  .  The  fact  that  internal  variables  of  the  lattice  are  entries  in  the  Cholesky 


factor  H  is  explained  in  3  .  We  present  the  generalizations  of  these  results  to  the  close 
to-Toeplitz  cases  in  the  next  section. 


VI.  Lattice  Presentation 


The  lattice  representation  of  the  coupled  recursions  in  the  correlation  method  of  linear 
prediction  yields  additional  insight.  The  transformation  from  reflection  coefficients  to 
correlation  values  is  particularly  easy  using  the  all-pole  (AR)  lattice  filter  7  .  With  a  zero 

input  sequence  and  an  initial  state  set  to  r<(,0 . 0  .  the  output  of  the  filter  is  the  causal 

part  of  the  correlation  sequence.  It  was  also  shown  in  3,  that  internal  variables  are  scaled 
entries  of  the  lower-triangular  matrix  H,  the  Cholesky  factor  of  the  Toeplitz  matrix  R. 
In  5i,  a  lattice  structure  of  more  then  two  lines  was  used  to  implement  the  normalized 
version  of  the  generalized  Leroux-Gueguen  algorithm,  for  all  the  cases  of  linear  prediction. 

Figure  3  showns  an  all-zero  (MA)  lattice  filter  that  may  be  used  to  implement  the 
recursions  of  the  previous  sections  for  the  covariance  case  of  linear  prediction.  The  other 
cases  correspond  to  simplified  versions  of  the  filter,  except  for  the  forward-backward  co 
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variance  case  which  requires  two  extra  lines.  The  conventions  used  in  Figure  3  are  the 
following  : 


K  A,  ’A'1 

*’(!)■ 

K'(2) 

K{  2) 

A'*(3) 

.  A' '  ( 3 ) . 

The  lattice  of  Figure  3  is  the  same  as  that  of  5  except  for  the  internal  organization  of  the 
cells  and  the  normalization  of  the  variables.  The  impulse  response  of  the  lattice  is  obtained 
by  placing  an  impulse  on  the  first  two  lines  and  zero  on  the  last  lines.  This  produces  the 
Generalized  Levinson  algorithm  of  section  III  (without  the  computation  of  the  reflection 
coefficients,  of  course).  In  Figure  4.  the  corresponding  all-pole  ( AR )  lattice  filter  is  shown. 
This  lattice  structure  has  one  very  interesting  application.  It  may  be  used  to  duplicate  the 
procedures  of  3  ,  7  .  8'  for  computing  correlations  in  the  close- to  Toeplitz  cases.  When  the 

inputs  is  zero  and  the  interna]  state  is  initialized  at  r(,  ,<».() . 0  .  the  first  row  (or  column) 

of  the  covariance  matrix  is  generated  on  the  first  output  line,  the  last  n  data  values  on 
the  second  line,  and  the  first  n  data  values  are  generated  on  the  third  line.  The  internal 
variables  at  the  inputs  (or  outputs)  on  the  first  line  of  the  cells  reproduce  the  entries  of 
the  Choleskv  factor  HD2.  The  entries  of  cell  ?  reproduce,  in  time  sequence,  the  entries  of 
the  column  of  H .  exactly  as  in  the  Toeplitz  case. 

Using  the  following  notation  for  the  entries  of  HD 2  : 


-H"(  0) 

tf°(l)  H'(  1)  <r> 

HD2  H" H'  •••//"  (17.2) 

Hn  ’(»  1) 

.H"(n)  Hl(n)  ...  A"  *(n )  Hn(n). 
then  the  inptit  of  cell  (?  -*  1)  (output  of  cell  i)  at  time  j  i  (on  the  top  line)  is 
Because  of  the  delay,  one  can  see  that  this  variable  will  be  zero  for  j  <  i.  Equivalently, 
entries  on  the  A-^1  line  are  equal  to  the  corresponding  entry  of  the  (A’  1)'^  column  of 

~  I 

H  .  This  is  illustrated  in  Figure  5.  The  algorithm  for  actually  computing  all  of  these 
variables  from  the  reflection  coefficients  is  a  generalization  of  that  given  in  7  .  We  may 
take  advantage  of  the  fact  that  almost  half  of  the  computed  variables  in  the  algorithm  of 
8  are  equal  to  zero,  as  Hl(j)  and  H  {j  -  1)  are  zero  for  j  <  i.  to  reduce  the  number  of 
computations.  The  algorithm  is  then 


Initialization  :  7/°(0)  -  rU  M 
For  j  1 .....  n  : 


<>r  i  j  1. 


H  (j)  -  H  (j)-H'(j  1 ) A',7!  , 

*Th\j) 


For  )  -  0. 


■J  1  : 


H'-'(j)  ^*0-  1)-  H'(j) 

If  the  algorithm  is  initialized  with  77u( 0)  -  1  instead  of  rUtU.  then  all  the  variables 
are  scaled  down  by  ru  o.  When  only  the  reflection  coefficients  are  known  then  both  the 
prediction  errors  ct *  and  A,  may  be  recovered  using  the  update  formula  given  in  section 
III.  Note  that  in  this  algorithm  H  (j)  is  the  row  of  H  .  which  is  a  row  vector  made  up 
of  (a  -  1)  elements.  Note  also  that  the  algorithm  is  not  restricted  to  the  methods  of  linear 
prediction  but  is  applicable  for  all  the  possible  values  of  the  displacement  rank  o.  When 
this  algorithm  is  run  in  conjunction  with  the  recursions  for  V'.  C  and  A',.  then  we  have  a 
complete  QR  algorithm  for  computing  V  and  H  in  the  QR  factorization  1  V H  1 . 

VII.  Multi-variable  case 

In  several  applications,  multi-variable  linear  prediction  is  necessary,  by  which  we  mean 
that  the  datum  y,  is  a  vector  of  length  d.  The  data  may  come  from  several  sensors  in  a 
multi-channel  system,  or  simply  be  a  collection  of  scalar  variables  in  one  row  (or  column) 
of  a  two  dimensional  image.  The  prediction  coefficients  a}  become  d  by  d  matrices.  Tin- 
forward  linear  prediction  problem  is  to  predict  the  value  of  y,  as 


N "  a}y,  j  for  ?  —  k . /. 


(17/.  1 


The  prediction  error  is 


-  y,  -  y,  ■-  y,  -  a,?/.  , 


(177.2) 


t 


X  .  A  .V  .N.V.VoS  ■jAj'V'A' 


which  can  still  be  organized  in  a  matrix  fashion  as  Y8 


(VII. 2) 


1  is  a  (/  A’  -  1)  by  (n  -t-  1  )d  block  Toeplitz  matrix.  Then,  depending  on  the  values  of  k 
and  /.  one  ends  up  with  different  methods  of  linear  prediction.  The  solution  of  the  problem 
is  the  same  as  in  the  scalar  case  except  for  element  dimensions  and  occasional  transpose. 
Our  purpose  in  this  section  is  just  to  give  the  generalization  of  our  vector  algorithms 
to  the  multi- %-ariable  case.  The  Grammian  matrix  R  is  now  block  close-to- Toeplitz  with 
ri,j  —  rJtl ■  The  generalized  Levinson  algorithm  performs  the  block  Cholesky  factorization 
of  R  1 ,  the  generalized  QR  algorithm  the  block  QR  factorization  of  Y .  and  the  generalized 
LeRoux-Gueguen  algorithm  the  block  Cholesky  factorization  of  R. 

The  generalized  Levinson  recursions  are 


B'+x  =  B'  -  Z  A'  Ay1  A',7;, 
A^1  =  Z  A ’  -i  B'  A/y  1  A'„  j 


for  i  -  0. 


n  -  1. 


(VII  A) 


with  the  initializations 


4°  =  7  0  ...  0 


ATo  —  rU(, 


[  c(). 0  0 

[  0  -s 


(VII. o) 


The  dimensions  of  the  variables  involved  are  the  following  :  .4'  is  (??  -  1  )d  by  d.  Bl  is 
(n  -+■  1  )d  by  (a  1  )d,  A',  is  d  by  d,  Mt  is  (a  —  1  )d  by  (o  -  1  )d.  S  is  ad  by  ad.  C  is  nd  by 
ad.  and  finally  Z  is  (n  -  l)<f  bv  (n  -t  1  )d. 

The  matrix  Z  in  these  equations  is  the  delay  matrix  : 


0  •  •  •  0  0 


(177.6) 


e.  ",  v.  .*  v,  „• .>  .*  a  .*  -■  .■  •  ,•  „•  ,«  v  •««  .<  v  v  v  'e  ,« •>  v  - j  v  ■.* 


A',.]  is  an  (a  -*  1  )d  by  d  matrix  and  is  computed  by 


A',7;,  =  -{ZA'Y 


(VII.  7 


The  prediction  error  matrix  X, .  and  the  error  matrix  A/,  are  updated  as  follows  : 

_  \T .  I \'T  t' 


AT-i  =  iV,  -  K,r  ,  A/,1  A',^1 
M-.  =  A/,  -  A\+,  N~'  Kl, 

The  lattice  recursions  for  the  matrix  columns  of  U  in  the  covariance  case  are 


(VII.  8) 


■“t-  1  ^-1 


5'  '  =  5  4  Z  S’  N-'hZ, 

S’* 1  =  Z  5’  -  5*  A/,  1  A',  ,  i  for  i  =  0, . . . ,  n  -  1. 


(r//.9) 


with  the  initializations 


yA-n-l 


JV«  -  r„,0  -  (t*0)7^0 


J/l,  0  0- 


,vL„  0  0 

yjl-n-  1  0  0 


y,f  o  o. 


A\,  o  0 

A7„  =070 
0  0  - 1 


(VI 1. 10) 


The  three  matrix  components  of  the  reflection  coefficient  A',,  j  are  computed  using 


*•£,< i)=-  i(£-")r.oijsi(.r) 

A’4  ,(2|  =  f"(-V  -  n }  ~  S'(.V) 
A'.r:  ,(3)  =  -S‘(n)  =  T» 


(1*77.12) 


The  other  forward  or  backward  cases  of  linear  prediction  are  simplified  version  of  this 
algorithm,  and  the  forward- backward  case  consists  of  two  of  these  recursions. 


WJSSWWS! SKW 


rvwii* 


%  .S  /»  *■  A  A  u\  -  -  .  •  A  .  *  m  *  A  -  »*  »  •  A 


1$ 

s 


I'Jt  »<«  .«> 


Finally,  the  Schur  recursions  are 

H'1  =  H  -  Z  H'  N~ 1  A,r  , 
H"1  =  Z  H'  4  tf'  M,  1  AV, 

with  the  initializations 


(!'//. 13) 


'  ^0,0  ■ 

r,,u 

-  u 

-  7"n,0  - 

H  = 

for  i  -  0. . .  . ,  n  1 . 


0  0 
n,0 


No  =  r0,  o 


ATn  0 


(V//.14) 


Note  that  Hl  is  the  block  column  of  HD2,  or  similarly  H'  Nt  1  is  the  block  column 
ofH.  D  is  the  block  diagonal  matrix  which  block  diagonal  element  is  Nt.  Note  also 

that  the  first  i  blocks  of  H'  and  the  first  (i  -  1)  block  rows  of  H  are  equal  to  zero.  The 
coefficient  A', +  j  may  be  read  out  of  the  recursions  as 


H  (i)  =  0 . o:  H  (t)  -  A'£, 


(V  11.1b) 


so  that  -  KJl ,  equals  the  first  non-zero  block  row  in  H  .  These  coupled  recursions  include 
the  update  for  N,  —  H'(i),  the  first  non-zer  >  element  in  H’.  The  only  necessary  side 
computation  is  the  update  of  M,  : 


M,.,  M,  A',-i  AV  1  A',7;, 


(17/.  16) 


Note  that  the  use  of  a  lattice  to  compute  the  Cholesky  factor  H  from  reflection  coefficients 
is  still  valid  in  the  multi-variable  case. 

VIII.  Conclusion 

We  have  derived  vector  algorithms  for  Cholesky  and  QR  factoring  Toeplitz  and  close- 
to-Toeplitz  matrices  for  all  of  the  cases  of  linear  prediction.  The  same  coupled  recursions 
are  used  in  all  the  algorithms,  namely 

V*41  =  N'  -  Z  M'K*, 

/l/141  =  Z  M'  4  V*  A"1  AY., 


Is? 


y 

T 

Jl! 

A 

tt 
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The  vector  M‘  contains  the  iin  column  of  the  matrix  A.  HD2.  or  U,  depending  upon  which 
factorization  is  being  computed.  The  inner  products  required  to  compute  the  reflection 
coefficients  and  to  initialize  the  variables  are 

•  Ml  ~  .4'  :  inner  product  required  for  computing  r,  „  and  K,. 

•  M'  =  V'  :  inner  product  required  for  computing  Kt  only. 

•  M 1  =  H'a 2  :  inner  product  required  for  computing  r,.0  only. 

The  Choleskv  algorithms  have  complexity  n2a.  The  fast  algorithms  for  the  orthogonal 
matrix  U  have  complexity  Nna,  where  N  is  the  number  of  data  values  available. 

The  QR  factorization  of  the  data  matrix  Y  in  the  covariance  case  of  linear  prediction 
may  be  used  for  a  general  Toeplitz  matrix  T.  This  gives  a  fast  algorithm  to  obtain  T  =  UR 
in  a  QR  method  to  obtain  eigenvalues.  If  Y  is  symmetric  (meaning  that  the  time  series  is 
symmetric  yn+t  =  yn-t)  then  the  forward  and  backward  predictors  are  reversals  of  each 
other,  or  in  other  words  U'  is  the  reversal  of  the  first  column  of  U  .  which  simplifies  the 
computation. 
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Appendix 

Fast  QR  algorithm  for  the  forward-backward  covariance  method  of  linear 
prediction. 

The  data  matrix  in  the  forward- backward  covariance  method  of  linear  prediction,  is 
the  following  : 


Yfb 


(.4.1) 


where  V  and  A*  are  the  forward  and  backward  data  matrices,  respectively,  in  the  covariance 
method  of  linear  prediction.  The  Grammian  is 


R  YjhYfh  =  YtY  +  _YrA*  -  YtY  4  JYTYJ 


(.4.2) 


which  means  that  R  is  a  centro-symmetric  matrix.  Its  displacement  rank  is  four,  and  the 


displacement  matrix  b'<R\  is  given  by 
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(.4.3) 


Note  that  f.R  -  -Jb  R  J.  As  the  generalized  Levinson  algorithm  does  not  use  the  centro- 
symmetry  property  of  R  (only  J2o  and  Rn  are  centro- symmetric),  the  matrix  .4  is  n  by  5. 
As  before,  the  recursions  for  A1  and  A  induce  recursions  for  the  orthogonal  matrix  in  the 
QR  factorization  : 


and  its  column 


1/6-4  = 


Yf„A'  = 
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A  = 


A‘  = 
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V * 
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(.4.4) 


(.4.5) 


Once  more  the  “generic”  recursions  have  to  be  used  to  get  recursions  for  the  columns  of 
the  orthogonal  matrix. 

Define  the  matrices  TV  and  >V  to  be  the  Toeplitz  extension  of  Y  and  X,  respectively  : 


W  = 


T 
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W  = 


T 
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(-4.6) 
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(.4. 


The  recursions  are  then  derived  exactlv  as  before  with 
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The  recursions  are 
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=  Z  |r)  +  A,  1  A',,  )  for  7  =  0, . . . ,  n  -  1 . 


with  the  initializations 
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The  five  components  of  the  reflection  coefficient  Al+i  are  computed  using 


<t,2AVi(1)  -  ,(lTV)r  ,0 


W[n) 

r* 


(r"  )T  .o 


V(n) 


<T,2A',4j(2)  =  r*(.V  -  n )  i  (3) 

fr,2AI+1(4)  =  -V(n)  cfA'i-.  ](5)  -  V*(.V  n) 


which  comes  from  the  use  of 


rU(l  0  0  0  0 

ri.u  -y0  yA-n  yn  i  -yA 


•  rn,u  -y,.  i  yA  i  yo  -y\ 


Thore  are  just  two  inner  products  per  update. 
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problems.  The  QR  factor  orthogonalizes  the  data  matrix  and  solves  the  problem  of  Cholesky 
factoring  the  experimental  correlation  matrix  and  its  inverse.  This  means  we  may  use  genera! 
Levinson  algorithms  to  derive  generalized  QR  algorithms,  which  are  then  used  to  derive  generi 
ized  Schur  algorithms.  All  three  algorithms  are  true  lattice  algorithms  that  may  be  implemen 
either  on  a  vector  machine  or  on  a  multi-tier  lattice,  and  all  three  algorithms  generate  gem 
ralized  reflection  coefficients  that  may  be  used  for  filtering  or  classification,  i 
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